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Abstract 

We  develop  an  alternative  formulation  of  conservative  finite  difference  weighted  essen¬ 
tially  non-oscillatory  (WENO)  schemes  to  solve  conservation  laws.  In  this  formulation, 
the  WENO  interpolation  of  the  solution  and  its  derivatives  are  used  to  directly  construct 
the  numerical  flux,  instead  of  the  usual  practice  of  reconstructing  the  flux  functions.  Even 
though  this  formulation  is  more  expensive  than  the  standard  formulation,  it  does  have 
several  advantages.  The  first  advantage  is  that  arbitrary  monotone  fluxes  can  be  used  in 
this  framework,  while  the  traditional  practice  of  reconstructing  flux  functions  can  be  ap¬ 
plied  only  to  smooth  flux  splitting.  The  second  advantage,  which  is  fully  explored  in  this 
paper,  is  that  it  is  more  straightforward  to  construct  a  Lax-Wendroff  time  discretization 
procedure,  with  a  Taylor  expansion  in  time  and  with  all  time  derivatives  replaced  by  spa¬ 
tial  derivatives  through  the  partial  differential  equations,  resulting  in  a  narrower  effective 
stencil  compared  with  previous  high  order  finite  difference  WENO  scheme  based  on  the 
reconstruction  of  flux  functions  with  a  Lax-Wendroff  time  discretization.  We  will  de¬ 
scribe  the  scheme  formulation  and  present  numerical  tests  for  one-  and  two-dimensional 
scalar  and  system  conservation  laws  demonstrating  the  designed  high  order  accuracy  and 
non-oscillatory  performance  of  the  schemes  constructed  in  this  paper. 
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1  Introduction 


In  this  paper,  we  are  interested  in  designing  an  alternative  formulation  of  high  order 
conservative  finite  difference  WENO  (weighted  essentially  non-oscillatory)  methods,  with 
the  Lax-Wendroff  time  discretizations,  in  solving  nonlinear  hyperbolic  conservation  laws 

f  ut  +  V  ■  f(u)  =  0  m 

\  u(x,  0)  =  u0{x). 

On  a  uniform  mesh  ay  =  iAx,  a  semi-discrete  conservative  finite  difference  scheme 
for  solving  (1)  has  the  following  form 


dt 


Ui 


1 

Ax 


fi+h 


(2) 


where  Ui  is  an  approximation  to  the  point  value  u(xi,t),  and  the  numerical  flux 


/«+  i  f  ( Ui—r ,  .  .  .  ,  Ui- |_s) 

is  designed  so  that 

^  (/i+i  -  /i-i)  =  f(u(x)) x\Xi  +  0(Axfc)  (3) 

for  a  fc-th  order  scheme.  Most  of  the  high  order  finite  difference  schemes,  for  example 
the  high  order  finite  difference  WENO  schemes  in  [4,  1,  12],  use  the  following  simple 
Lemma  in  [14]: 


Lemma  1.1  (Shu  and  Osher  [14]):  If  a  function  h(x)  satisfies  the  following  relationship 


/ OOe>)  =  -r-  /  (4) 

Ax  Jx_  A, 

then 

1  f  /\t  At  \ 

^  +  -  2")  -  h(x  -  — )J  =  /(«(*))..  (5) 

The  proof  of  this  lemma  is  a  simple  application  of  the  fundamental  theorem  of  cal¬ 
culus. 
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With  this  simple  lemma,  it  is  clear  that  we  should  take  the  numerical  flux  as 

fi+i  =  h(xi+i)  +  0(Axk)  (6) 

with  a  smooth  coefficient  in  the  0( Axk)  term,  to  obtain  the  high  order  accuracy  (3). 
Since  the  definition  of  the  function  h(x)  in  (4)  implies 

1  fx‘+i 

hi  =  -^J  =  /te) 

which  is  known  once  we  have  the  point  values  iti,  we  are  facing  the  standard  recon¬ 
struction  problem  of  knowing  the  cell  averages  {hi}  of  the  function  h(x),  and  seeking  to 
reconstruct  its  point  values  at  the  cell  boundaries  h(xi+i)  to  high  order  accuracy,  to  be 
used  as  the  numerical  flux  fi+ i  in  (6).  This  reconstruction  problem  is  the  backbone  of 
the  standard  finite  volume  schemes,  for  example  the  essentially  non-oscillatory  (ENO) 
finite  volume  schemes  [3]  and  WENO  finite  volume  schemes  [7].  Therefore,  using  the 
simple  Lemma  1.1,  we  can  simply  apply  the  standard  reconstruction  subroutine  in  a 
high  order  finite  volume  scheme  to  obtain  the  numerical  flux  {/,:+i }  for  a  high  order 
conservative  finite  difference  scheme  from  the  flux  point  values  {/(«,)}. 

Because  of  the  clean  conception  and  easy  implementation,  this  formulation  of  numer¬ 
ical  fluxes  in  high  order  conservative  finite  difference  schemes  has  been  used  widely,  for 
example  in  high  order  finite  difference  ENO  schemes  [14]  and  WENO  schemes  [4,  1,  12]. 
However,  there  are  several  disadvantages  of  this  formulation: 

1.  In  order  to  achieve  upwinding  and  nonlinear  stability,  an  approximate  Riemann 
solver  or  a  two-point  monotone  flux  is  used  in  finite  volume  schemes.  With  this 
finite  difference  flux  formulation,  upwinding  can  be  implemented  easily  only  for 
the  following  smooth  flux  splitting 

f(u)  =  /+0)  +  f~(u) 

where  -£;f+(u)  >  0  (or  in  system  case,  this  Jacobian  matrix  has  only  real  and 
positive  eigenvalues)  and  <  0  (or  in  system  case,  this  Jacobian  matrix  has 
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only  real  and  negative  eigenvalues).  In  order  to  guarantee  high  order  accuracy,  the 
two  functions  f+(u )  and  f~(u )  should  be  as  smooth  functions  of  u  as  f{u).  The 
most  common  flux  splitting  used  in  finite  difference  schemes  is  therefore  the  Lax- 
Friedrichs  flux  splitting,  which  is  the  most  diffusive  among  two-point  monotone 
fluxes. 

2.  The  Lax-Wendroff  time  discretization  procedure  [8]  is  difficult  to  design  and  results 
in  a  rather  wide  effective  stencil.  We  will  elaborate  in  more  details  on  this  point 
later  in  the  paper. 

3.  Since  the  reconstruction  is  performed  directly  on  the  flux  values  {/(tq)}  (or  {f+(ui)} 
and  {f~(ui)}  with  flux  splitting),  not  on  the  point  values  of  the  solution  {ui},  it 
is  more  difficult  to  maintain  free  stream  solutions  exactly  in  curvilinear  meshes 
for  multi-dimensional  flow  computation.  This  is  because  the  fluxes  in  curvilinear 
coordinates  involve  metric  derivatives,  resulting  in  non-exact  cancellations  when 
nonlinear  reconstructions  are  performed  for  different  fluxes. 

We  would  like  to  explore  an  alternative  formulation  for  constructing  numerical  fluxes 
in  high  order  conservative  finite  difference  schemes,  originally  designed  in  [13],  which 
involves  interpolations  directly  on  the  point  values  of  the  solution  {iq}  rather  than  on 
the  flux  values.  This  alternative  formulation,  even  though  less  clean  and  more  computa¬ 
tionally  expensive,  does  overcome  all  three  disadvantages  listed  above.  We  will  explore 
the  overcoming  of  the  first  two  disadvantages,  especially  the  second  one,  in  detail  in  this 
paper. 

We  now  give  a  brief  survey  of  WENO  schemes  and  Lax-Wendroff  time  discretiza¬ 
tion.  WENO  schemes  are  high  order  schemes  for  approximating  hyperbolic  conservation 
laws  and  other  convection  dominated  partial  differential  equations.  They  can  produce 
sharp,  non-oscillatory  discontinuity  transitions  and  high  order  accurate  resolutions  for 
the  smooth  part  of  the  solution.  The  first  WENO  scheme  was  introduced  in  1994  by 
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Liu,  Osher  and  Chan  in  their  pioneering  paper  [7],  in  which  a  third  order  accurate  fi¬ 
nite  volume  WENO  scheme  in  one  space  dimension  was  constructed.  In  [4],  a  general 
framework  is  provided  to  design  arbitrary  high  order  accurate  finite  difference  WENO 
schemes,  which  are  more  efficient  for  multi-dimension  calculations.  Very  high  order  hnite 
difference  WENO  schemes  are  documented  in  [1]. 

The  Lax-Wendroff  time  discretization  procedure,  which  is  also  called  the  Taylor  type 
time  discretization,  is  based  on  the  idea  of  the  classical  Lax-Wendroff  scheme  [5],  relying 
on  writing  the  solution  at  the  next  time  step  by  Taylor  expansion  in  time,  and  repeatedly 
using  the  partial  differential  equation  (PDE)  to  rewrite  the  time  derivatives  in  this  Taylor 
expansion  as  spatial  derivatives.  These  spatial  derivatives  can  then  be  discretized  by 
standard  approximation  procedures,  for  example  the  WENO  approximation.  We  denote 
■u™  as  the  approximation  of  the  point  values  u(xi,tn),  and  u^r>  as  the  r-th  order  time 
derivative  of  u,  namely  =  ^f.  We  also  use  u u"  and  u!"  to  denote  the  first  three 
time  derivatives  of  u.  By  Taylor  expansion,  we  obtain 


u(x ,  t  +  At)  =  u(x ,  t )  +  Atu'  + 


At2 

i>r 


u"  + 


At 3  At4 

■u  + 


3! 


4! 


«<4>  +  + 

5! 


(7) 


If  we  approximate  the  first  k  derivatives,  we  can  get  k- th  order  accuracy  in  time.  The 
main  idea  of  the  Lax-Wendroff  time  discretization  procedure  is  to  rewrite  the  time  deriva¬ 
tives  vx’  as  spatial  derivatives  utilizing  the  PDE,  then  discretize  these  spatial  derivatives 
to  sufficient  accuracy.  This  is  rather  straightforward  for  a  hnite  volume  scheme.  However, 
for  a  hnite  difference  scheme  using  the  formulation  (6)  and  Lemma  1.1  to  compute  its 
numerical  huxes,  it  is  more  difficult  to  implement  the  Lax-Wendroff  procedure  because 
we  do  not  have  approximations  to  the  spatial  derivatives  of  the  solution  u  at  our  disposal, 
only  approximations  to  the  spatial  derivatives  of  the  hux  function  f(u).  In  [8],  Qiu  and 
Shu  designed  a  hnite  difference  WENO  scheme  with  Lax-Wendroff  time  discretization 
using  the  formulation  (6)  and  Lemma  1.1.  To  avoid  the  difficulty  mentioned  above,  they 
use  the  approximations  to  the  lower  order  time  derivatives  to  approximate  the  higher 
order  ones.  As  a  consequence,  the  method  in  [8]  involves  a  rather  wide  effective  stencil. 


5 


We  will  provide  more  details  to  illustrate  this  phenomenon.  In  the  procedure  to  build  a 
fifth  order  in  space  and  fourth  order  in  time  scheme  in  [8]  for  solving  the  one-dimensional 
scalar  equation  ut  +  f(u)x  =  0,  the  following  steps  are  involved: 

Step  1.  Using  the  PDE,  we  obtain  v!  =  —f(u)x,  and  the  first  order  time  derivative 
can  be  obtained  by  the  conservative  finite  difference  WENO  approximation  (3).  To 
obtain  a  fifth  order  WENO  approximation  for  u' |™,  we  would  need  the  point  values 
S  =  {u]_ 3,---  ,u]+3}. 

Step  2.  When  constructing  the  second  order  time  derivative  u"  =  —{f'[u)u')Xl  we  denote 
()j  =  f'(uj)u'j,  where  u!-  are  the  point  values  of  the  first  order  time  derivative  computed 
in  the  previous  Step  1.  Then  we  can  use  a  fourth  order  central  difference  formula  to 
approximate  u"  =  —gx  at  the  point  (xj,tn): 

u"  =  “  8^— i  +  8ft'+ 1  -  9j+ 2)- 

It  was  observed  in  [8]  that  the  more  costly  WENO  approximation  is  not  necessary  for 
second  and  higher  order  time  derivatives,  as  the  standard  central  difference  approxima¬ 
tions  perform  well.  Considering  gj  =  f'(uj)u'j  and  using  Step  1  to  obtain  w',  we  have 
the  effective  stencil  for  obtaining  the  point  value  u" |”  to  be  S  =  {u"_5,  ■  •  •  ,w"+5}. 

Step  3.  When  constructing  the  third  order  time  derivative  u’"  =  —  (f,(u)u"+f'(u)(u')2)x, 
we  let  gj  =  ( f\uj)u"  +  /,/(uJ)(u,)|),  where  the  point  values  v!:]  and  u”  are  obtained  in 
the  previous  Step  1  and  Step  2,  respectively.  Then  we  repeat  Step  2  to  approximate  the 
third  order  time  derivative  using  a  fourth  order  central  difference  approximation  (third 
order  accuracy  is  enough  here,  but  central  approximations  provide  only  even  orders  of 
accuracy).  Considering  u'  and  u"  are  obtained  in  Step  1  and  Step  2,  we  obtain  the 
effective  stencil  for  approximating  v!"  at  the  point  ( Xj,tn )  to  be  S  =  {«"_ 7,  •  •  •  ,Wj+7}. 

Step  4.  The  construction  of  the  fourth  order  time  derivative  u ^  =  —(f(u)u'"  + 
3/"(!i)iiVT/",(«)(M,)3)i  is  obtained  in  a  similar  fashion.  Let  gj  =  f'(uj)uj,+3f"(uj)u'ju''+ 
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where  the  point  values  u'p  u"  and  u"-  are  obtained  in  the  previous  Steps 
1,  2  and  3  respectively.  Then  we  use  a  second  order  central  difference  to  approximate 
uf]  =  ~  2Kx(9i+1  ~  Thus,  to  get  the  value  at  the  point  ( Xptn ),  we  actually 

need  point  values  in  the  effective  stencil  S  =  {m™_8,  •  •  •  ,  8}. 

In  summary,  to  obtain  a  fifth  order  in  space  and  fourth  order  in  time  Lax-Wendroff 
type  WENO  scheme  in  [8],  the  effective  stencil  is  S  —  {/,_8,  •  ■  ■  ,  /j+8}  for  the  point  Xj, 
containing  17  grid  points.  This  is  a  pretty  wide  stencil.  Of  course,  it  is  still  narrower 
than  Runge-Kutta  time  discretization.  When  a  four-stage,  fourth  order  Runge-Kutta 
time  discretization  is  used  on  a  fifth  order  finite  difference  WENO  scheme,  the  width  of 
the  effective  stencil  is  25  grid  points. 

One  major  motivation  to  discuss  the  alternative  formulation  of  high  order  conser¬ 
vative  finite  difference  schemes  in  this  paper  is  to  make  the  effective  stencil  narrower 
when  applying  the  Lax-Wendroff  time  discretization.  We  will  review  WENO  interpola¬ 
tion  in  Section  2,  which  is  the  main  approximation  tool  for  the  alternative  formulation 
of  high  order  conservative  finite  difference  schemes.  Details  on  the  construction  and 
implementation  of  the  alternative  formulation  will  be  described  in  Section  3,  for  one- 
and  two-dimensional  scalar  and  system  equations.  In  Section  4,  extensive  numerical 
examples  are  provided  to  demonstrate  the  accuracy  and  effectiveness  of  the  methods. 
Concluding  remarks  are  given  in  Section  5. 

2  WENO  interpolation 

A  major  building  block  of  the  alternative  formulation  of  high  order  conservative  finite 
difference  scheme  discussed  in  this  paper  is  the  following  WENO  interpolation  procedure. 
Given  the  point  values  tp  =  u(xj)  of  a  piecewise  smooth  function  u(x),  we  would  like 
to  find  an  approximation  of  u(x)  at  the  half  nodes  xi+i.  This  WENO  interpolation 
procedure  has  been  developed  in,  e.g.  [10,  2,  12]  and  will  be  described  briefly  below. 

The  first  component  of  the  general  WENO  produce  is  a  set  of  lower  order  approxi- 
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mations  to  the  target  value,  based  on  different  stencils,  which  are  referred  to  as  the  small 
stencils.  We  would  find  a  unique  polynomial  of  degree  at  most  k  —  1,  denoted  by  pr(x), 
which  interpolates  the  function  u(x),  namely  pr{.Xj)  =  Uj,  at  the  mesh  points  Xj  in  the 
small  stencil  Sr  =  {xi-r,  Xj_r+i, . . . ,  Xi+S }  for  r  =  0,1 , ,k  —  1,  with  r  +  s  +  1  =  k. 
Then,  we  would  use  =  pr(xi+i)  as  an  approximation  to  the  value  u(xi+ i),  and  we 
have 

«•+ 1  -«(xi+i)  =  0(Axfc)  (8) 

if  the  function  u(x)  is  smooth  in  the  stencil  Sr.  For  example,  when  k  =  3,  we  have 

SO  li^h+2}:  Si  \^Xi—  i,  Xi:  S2  2,  %i— 1,  aild 


(0)  3  3  1 

=  +  Tl+1  -  8Ui+2’ 


1  3  3 

—  gM*-i  S  +  -Wj+i, 


(2)  3  5  15 

•A  =  s“'-2 "  J”*-1  +  T“- 


The  second  component  of  the  general  WENO  procedure  is  to  combine  the  set  of  lower 
order  approximations  to  form  a  higher  order  approximation  to  the  same  target  value, 
based  on  the  big  stencil  S  =  {ay_fc+ 1,  •  •  •  ,£*+&- 1},  which  is  the  union  of  all  the  small 
stencils  in  the  first  component.  We  would  find  a  unique  polynomial  of  degree  at  most 
2k  —  2,  denoted  by  p(x),  interpolating  the  function  u{x)  at  the  points  in  the  big  stencil. 
We  use  ui+i  =  p(xi+ 1)  as  the  approximation  to  the  value  u(xi+ 1),  and  we  have 


m(+i  -n(xt+i)  —  0(Ax2k  ‘) 


(9) 


provided  that  the  function  u{x)  is  smooth  in  the  big  stencil  S.  When  k  =  3,  we  have 


3 

2  128 ' 


5 

32 


^i+1  2  t  r.  A  ^2  ~b  00  ^4- 2* 


45  15 

— Ui  H - i 

64  32 


5 

128 


There  are  constants  dr  such  that 


k- 1 


uiH 


drU\+b 


r= 0 


(10) 


where  the  constants  satisfy  Ylr= o  (k  =  1-  These  constants  are  usually  referred  as  the 
linear  weights  in  the  WENO  procedure.  For  the  case  k  —  3,  the  linear  weights  are  given 
as 


do  — 


5 

16’ 


G?2  — 


1 

16 


The  third  component  of  the  general  WENO  procedure  is  to  change  the  linear  weights 
into  nonlinear  weights,  to  ensure  both  accuracy  in  smooth  cases  and  non-oscillatory 
performance  when  at  least  one  of  the  small  stencils  contain  a  discontinuity  of  the  function 
u(x).  We  can  achieve  this  through  changing  the  linear  weights  dr  into  the  nonlinear 
weights  ay,  defined  by 


ar 


ay  = 


\~~\k—l 
a^s=o  a‘ 


-,  r  —  0, . . . ,  k  —  1 


(11) 


with 


dr 


ay  = 


(12) 


(e  +  (3r)2 

Here,  we  take  e  =  10-6,  avoiding  the  denominator  to  be  zero,  and  (3r  are  the  so-called 
smoothness  indicators  of  the  stencil  Sr,  which  measures  the  smoothness  of  u(x)  in  this 
stencil.  Following  the  original  WENO  recipe  [4],  we  define 


k  1  rt-j* 

1=1  Jxi-  A 


.  21—1  ( dlPr(x) 


Ax 


dlx 


dx. 


When  k  =  3,  the  formula  gives 


13  1 

Po  —  ~(M*  —  ^ Ui+ 1  T  Ui+2Y  +  ^(3 Ui  —  4rtj_)_i  +  Ui+2)2 , 

13  1 

Pi  —  Y2yUi~1  ~  ^ Ui  Ui+ 1)2  T  1  —  ui+l)2 1 

13  1 

P2  —  —  2rtj_i  +  Ui)2  +  —  (rtj_  2  —  4  Ui~i  +  3rtj)2. 

Thus,  we  get  the  WENO  interpolation  of  u(x)  at  the  point  xi+i  as 


(13) 


fc-i 


ui+b  = 


E(r) 

^Ui+  i- 


r=0 


(14) 
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We  denote  it  by  u  x  since  the  big  stencil  S  used  to  obtain  this  approximation  is  biased 
*+2 

to  the  left  (there  is  one  more  grid  point  to  the  left  of  xi+i  than  to  the  right).  If  u(x)  is 
smooth  in  this  big  stencil  S,  we  have 

-  u(xi+ 1)  =  0( Ax2*-1). 

If  one  of  the  small  stencils  contains  a  discontinuity,  the  corresponding  smoothness  indica¬ 
tor  will  be  much  larger  than  that  of  other  small  stencils,  hence  the  nonlinear  weight  corre¬ 
sponding  to  this  small  stencil  will  be  very  small,  thus  achieving  essentially  non-oscillatory 
performance.  If  we  choose  the  big  stencil  as  S'  =  {xi-k+ 2,  •  •  •  ,Xi+k}  which  is  biased  to 
the  right,  and  the  small  stencils  as  Sr  =  {xi-r+ 1,  ■  •  • ,  £;-r+fc},  for  r  =  0, . . . ,  k  —  1,  we 
can  obtain  in  a  similar  fashion  a  WENO  interpolation  of  u(x)  at  the  point  xi+i,  denoted 

by  u+  1,  which  has  the  same  high  order  accuracy  in  smooth  regions  and  essentially  non- 
*+2 

oscillatory  performance  near  discontinuities  as  u7+ I11  fact,  the  process  to  obtain  u++ 1 
is  mirror-symmetric  to  that  for  u7+ 1,  with  respect  to  the  target  point  xi+i. 

3  Construction  of  the  scheme 

3.1  Construction  of  the  numerical  flux  in  one-dimension 

Assuming  that  f(u)  is  a  smooth  function  of  u,  we  would  like  to  find  a  consistent  numerical 
flux  function 

fi+±=f{ui—ri---iui+s)i  (15) 

such  that  the  flux  difference  approximates  the  derivative  f(u(x))x  to  k- th  order  accuracy 

-  fi-i)  =  f(u{x))x\xi  +  0{ Axk).  (16) 

The  numerical  flux  should  satisfy  the  following  conditions: 

•  /  is  a  Lipschitz  continuous  function  in  all  the  arguments; 

•  /  is  consistent  with  the  physical  flux  /,  that  is,  f(u,  •••,«)  =  f{u). 
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It  has  been  shown  in  [13]  that  there  exist  constants  a2,  a4, . . .  such  that 


Li  =/<+J+  V  a2kAx2k  (  )  +0{ AV«)  (17) 

fc=l  '  '  *+3 

which  guarantees  k  =  r-th  order  accuracy  in  (16).  The  coefficients  a2r  in  (17)  can 
be  obtained  through  Taylor  expansion  and  the  accuracy  constraint  (16).  To  get  an 
approximation  with  fifth  order  accuracy  {k  —  5  in  (16)),  we  can  use  the  first  three  terms 
given  as 

-  1  7 

h+ \  =  fi+ 1  —  24^X  5700 fxxxx\i+h'  (l®) 

The  first  term  of  the  numerical  flux  in  (18)  is  approximated  by 


fi+i  = 


(19) 


with  the  values  uf  x  obtained  by  the  WENO  interpolation  discussed  in  the  previous 
l+2 

section.  The  two-argument  function  h  is  a  monotone  flux.  It  satisfies: 


•  h(a,  b)  is  a  Lipschitz  continuous  function  in  both  argument; 

•  h(a,  b )  is  a  nondecreasing  function  in  a  and  a  nonincreasing  function  in  b.  Symbol¬ 
ically  h(T, |); 

•  h(a,b )  is  consistent  with  the  physical  flux  /,  that  is,  h(a,a)  =  /(a). 


Examples  of  the  monotone  fluxes  include: 


1.  The  Godunov  flux: 


f  mina<x<6/(it),  if  a  <b 
\  max6<x<J(«),  if  a>b 


(20) 


2.  The  Engquist-Osher  flux: 


(u),0)du  +  /  min(//(u),  0)du  +  /(0)  (21) 

Jo 


r 

h(a ,  b)  —  max(/' 

Jo 
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3.  The  Lax-Friedrichs  flux: 


h{a,b)  =  ^[/(a)  +  f(b)  -a(6-a)]  (22) 

where  a  =  maxu  |/'(u)|  is  a  constant.  The  maximum  is  taken  over  the  relevant 
range  of  u. 

More  discussions  on  monotone  fluxes  can  be  found  in,  e.g.  [6]. 

The  remaining  terms  of  the  numerical  flux  in  (18)  have  at  least  Ax2  in  their  co¬ 
efficients,  hence  they  only  need  lower  order  approximations  and  they  are  expected  to 
contribute  much  less  to  spurious  oscillations.  Therefore,  following  the  practice  in  [8], 
we  approximate  these  remaining  terms  by  simple  central  approximation  or  one-point 
upwind-biased  approximation  with  suitable  orders  of  accuracy,  without  using  the  more 
expensive  WENO  procedure. 

3.2  One-dimensional  scalar  case 

Consider  the  one-dimension  scalar  conservation  law 

/  ut  +  f(u)x  =  0,  ,  , 

\  u(x,0)=uo(x).  [Z6) 

We  describe  the  Lax-Wendroff  time  discretization  to  get  k- th  order  accuracy  in  time.  In 
this  paper,  we  will  implement  the  procedure  up  to  fourth  order  in  time. 

1.  The  approximation  of  the  first  time  derivative  u'  =  —f(u)x  follows  the  standard 
WENO  procedure.  We  construct  the  numerical  flux  with  (2r  — l)-th  order  accuracy 
as  described  in  the  previous  section.  In  particular,  the  points  {xi_r+ 1, . . .  ,xi+r} 
are  required  to  obtain  the  WENO  interpolation  uf+i-  For  the  spatial  derivatives 
at  xi+i,  we  use  (2r  —  1  —  fc)-th  order  simple  finite  difference  approximations 
based  on  the  same  points  {xj_r+i, . . . ,  Xj+r}.  In  our  numerical  experiments,  we  aim 
for  fifth  order  spatial  accuracy  and  take  r  =  3.  The  numerical  flux  corresponding 
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to  u'  is  then  given  as 


x2(f"(u)(ux)2  +  f'(u)uxx)\x.+  ,  +  x*(f""(u)(ux)4 

6  f  ('u)('U3;)  UXx  T  4  f  (u)uxUxxx  T  3/  ( u)(uxx )  T  f  (u)Uxxxx)  lx.  , 


such  that 

i)  _  /(i) 

Jj+i  I 

2Ax  2  =-^U  +  0(Ax4) 

where  i  in  the  h(uT  i,w+  i)  term  are  obtained  through  the  fifth  order  WENO 
interpolation,  and  the  values  of  u,  ux,  uxx,  uxxx  and  uxxxx  at  xi+i  in  the  remaining 
terms  are  obtained  from  the  same  polynomial  interpolation  based  on  the  points 
{Xi-2, • •  • ,  xi+3}. 


2.  The  next  step  is  to  obtain  the  approximation  of  the  second  time  derivative  u"  = 
—  (/'(«)«%  =  (( f{u))2ux)x .  We  let  g  =  —(f'(u))2ux,  and  repeat  the  process  in 
Step  1.  The  two  major  differences  from  Step  1  are:  (i)  the  approximation  can 
be  one  order  lower  in  accuracy  since  there  is  an  extra  At  in  the  coefficient;  and 
(ii)  there  is  no  need  to  involve  WENO  in  the  approximation,  also  because  of  the 
extra  At  in  the  coefficient.  Therefore,  we  would  still  only  need  the  interpolation 
polynomial  based  on  the  points  {xj_2, . . .  ,  Xj+3}  for  all  the  values  of  u,  uxi  uxx, 
uxxx  and  uxxxx  at  xi+i  in  the  following  expression 


f(2)  —  a\  -  ———a  I 
j  a\  1  y  \x . .  i  0  A  yxx \x . .  i 

z+2  24  l+2 

A  7*2 

=  ~{f\u)fux\Xi+h  +  —(2  U"{u)f{uxf  +  2  r(u)f\u){uxf 
+  6  f"{u)f\u)uxuxx  +  ( f(u))2uxxx)\x 


such  that 


m 


Ax 


=  -u”\Xi  +  0(Aa:4). 
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3.  The  approximation  of  the  third  time  derivative  v!"  =  — (3/"(n)(/,(n))2n^ 
+  {f\u)fuxx)x  is  obtained  in  a  similar  fashion.  Let  g  =  (3 
ul  +  if{u))3uxx),  we  can  obtain  the  numerical  flux  as 


Ax2 

=9\xi+i  ~  ~2^9xx \x.+  i 

Ar2 

=(3  f"{u){f\u)Y(uxy  +  u\u)fuxx)\Xi+h  -  —(3.  r\u){f\u))\uxy 

+  18/>)/''(u)/'(u)(«*)4  +  18  f"\u){f\u))\uxfuxx  +  6(/"  («))3(^)4 
+  36  (f"(u))2f(u)(ux)2uxx  +  12  f"{u){f{u)fuxuxxx  +  9  f\u){f{u)f{uxx)2 
T  (,/  (m))  Uihi)|i  i 

*+TI 


where  the  values  of  u,  ux,  uxx ,  nxxx  and  uxxxx  at  xi+i  can  all  be  obtained  from  the 
same  interpolation  polynomial  based  on  the  points  {xj_2, . . .  ,xi+3}.  We  will  then 
have 


m 

Jij-i 


£  (3) 
I 

1  2 


Ax 


=  -uw|Xi  +  0(Ax3) 


4.  The  final  step  is  the  approximation  of  the  fourth  time  derivative  w*-4’  =  (4  f'"(u)(f'(u)) 


{uxf  +  12 (f"(u))2(f'(u))2(ux)3  +  12 f'\u)(f\u)fuxuxxx  +  {f'{u))Auxxx)x.  Let  g  = 

-(4rX«)(/X«))3K)3+12(/"(n))2(f(«))2K)3+12/"(w)(//(«))3^xo:+(//(«))4fe), 


we  can  obtain  the  numerical  flux  as 


=  -  (4 f"\u){f\u))\uxf  +  12 U'\u))\f\u))\uxf  +  12 f"{u){f{u)fuxu: 

+  (f{u))*Uxxx)\x..i 


where  the  values  of  u,  ux ,  uxx ,  nxxx  and  uxxxx  at  xi+i  can  all  be  obtained  from  the 
same  interpolation  polynomial  based  on  the  points  {ay_ 2, . . .  ,  x;+3}.  We  will  then 
have 


f(4)  _  ^(4) 

Ax 


=  —  it^l 


+  0(Ax2 
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5.  We  can  form  the  numerical  flux  as  the  sum  of  all  the  fluxes  above 


and  obtain  the  final  scheme 


which  is  fifth  order  accurate  in  space  and  fourth  order  accurate  in  time. 

Of  course,  we  can  also  increase  the  time  accuracy  to  fifth  order,  however  heavier 
algebra  will  be  involved.  In  practice,  the  time  accuracy  can  usually  be  lower  than  spatial 
accuracy  to  achieve  the  desired  resolution  on  modest  meshes.  This  is  also  verified  in  our 
numerical  results  in  Section  4. 

Clearly,  the  fifth  order  in  space,  fourth  order  in  time  Lax-Wendroff  finite  difference 
scheme  constructed  above  has  an  effective  stencil  of  only  7  points,  the  same  as  the  semi¬ 
discrete  scheme.  This  compares  quite  favorably  with  the  Lax-Wendroff  type  scheme 
in  [8]  which  has  an  effective  stencil  of  17  points,  or  a  fourth  order  Runge-Kutta  finite 
difference  WENO  scheme  which  has  an  effective  stencil  of  25  points. 

3.3  One-dimensional  systems 

For  a  one- dimensional  system  of  conservation  laws  (23), 


u(x,  t )  =  (n1(x,  t),  ■  ■  ■  ,  um(x,  t))T 


(24) 


is  a  vector  and 


(25) 


is  a  vector  function  of  u.  Suppose  the  system  is  hyperbolic,  i.e.  the  Jacobian  f'(u )  has 
m  real  eigenvalues 


Ai(u)  <  . . .  <  A m{u) 


(26) 
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and  a  complete  set  of  independent  eigenvectors 


ri(u),...,rm(u). 


(27) 


We  denote  the  matrix  whose  columns  are  eigenvectors  by 


R(u)  =  (ri(u), . . . , rm(u)). 


(28) 


Then 


R  1(u)f'{u)R{u)  =  A  (u) 


(29) 


where  A (u)  is  the  diagonal  matrix  with  Ai(w), . . . ,  Am(w)  on  the  diagonal. 

Similar  to  the  scale  case,  we  replace  the  time  derivatives  by  spatial  derivatives  using 
the  PDE.  We  still  use  the  numerical  flux  for  each  component 


k= i 


when  approaching  gx. 

When  constructing  the  first  part  of  numerical  flux  for  ut,  i.e.  h(u7  1}ut  x)  in  / fy, , 


we  could  use  the  WENO  interpolation  procedure  on  each  component  of  u  as  in  the 
scalar  case.  However,  for  strong  shocks  this  procedure  may  lead  to  oscillatory  results. 
Interpolation  in  the  local  characteristic  fields,  even  though  more  costly,  may  effectively 
eliminate  such  spurious  oscillations.  The  procedure  of  interpolation  in  the  local  charac¬ 
teristic  fields  can  be  summarized  as  follows.  At  each  point  xi+i  where  the  computation 
of  the  numerical  flux  is  needed,  we  perform  the  following  steps. 

1.  Compute  an  average  state  ui+ 1,  using  either  the  simple  arithmetic  mean 


ui+ 1  —  g  (Ui  +  m»+i)j 


or  a  Roe  average  [9]  satisfying 


f{ui+ i)  -  f(ui)  =  f(ui+i)(ui+1  -  m) 


if  it  is  available. 


16 


2.  Compute  the  right  eigenvectors  and  the  left  eigenvectors  of  the  Jacobian  f'{ui+ 1), 
and  denote  them  by 

R  =  R{ui+ 1),  R”1  =  R~l(ui+ 1). 

3.  Transform  all  the  point  values  which  are  in  the  potential  stencil  of  the  WENO 
interpolation  for  the  obtaining  u±+1,  to  the  local  characteristic  fields.  For  our  fifth 
order  scheme,  we  need  to  compute 

Vj  —  R~lUj,  j  =  i  —  2, i  +  3.  (31) 


4.  Perform  the  scalar  WENO  interpolation  procedure  for  each  component  of  the  char¬ 
acteristic  variable  v,  to  obtain  the  corresponding  component  of  the  interpolation 


5.  Transform  back  into  physical  space 


For  the  two-argument  flux  h(u~,u+),  instead  of  the  monotone  fluxes  for  the  scalar 
case,  we  can  choose  from  a  variety  of  numerical  fluxes  based  on  exact  or  approximate 
Riemann  solvers  [15].  For  example,  we  can  use  the  Godunov  flux,  the  Lax-Friedrichs 
flux,  the  HLLC  flux,  etc.  This  is  one  of  the  advantages  of  the  alternative  formulation 
discussed  in  this  paper. 

For  all  remaining  terms  in  the  Lax-Wendroff  process,  we  can  use  simple  interpola¬ 
tion  polynomials  based  on  the  same  stencil  {xj_ 2, . . .  ,£,+3}  component-wise  to  obtain 
approximations  to  the  values  of  u,  ux,  uxx ,  uxxx  and  uxxxx  at  xi+i.  No  characteristic 
decomposition  or  WENO  is  needed.  We  remark  that  in  the  algebraic  manipulations  of 
the  Lax-Wendroff  process,  the  derivatives  f'(u),  f"(u),  f”'(u )  etc.  will  become  matrices 
and  higher  order  tensors,  increasing  the  algebraic  complexity  of  the  final  formulas. 
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3.4  Two-dimensional  cases 


Consider  the  two-dimensional  conservation  laws: 

f  ut  +  f(u)x  +  g(u)y  =  0, 

1  u(x,y,  0)  =  uq(x,  y). 

The  same  Taylor  expansion  still  gives 

u(x,  y,t  +  At)  =  u(x,  y,  t )  +  Atu  +  — - u  +  — —u  +  -——u t4)  H - 

2  6  24 


(32) 


(33) 


and  we  would  still  like  to  use  the  PDE  to  replace  time  derivatives  by  spatial  derivatives. 
The  first  time  derivative  v!  =  —f(u)x  —  g(u)y  can  be  approximated  in  a  dimension-by¬ 
dimension  fashion  as  in  standard  finite  difference  schemes  [13] .  For  example,  for  the  fifth 
order  case,  we  have 


Ax 2 


7  Ax4 


fi+ij  =hi(ui^ut+lj)  -  -^(/"MK)2  +  f\u)uxx) \(x.+i,yj)  +  -(f""(u)(uxy 

Qf  U xx  3  f  {Vj){VjX x)  4/*  {^j^x^xxx  f  (y^'U'xxxx)  \  (x .  i  ,yj) 

=h^\j+^Utj+^  ~  ^-(9"(u)(uy)2  +  g'(u)uyy)\(Xi,y.+  i)  +  g""(u)(uy )4 

+  6(7  (u)(uy)  Uyy  ?>g  (u)  (Uyy)  ~\~4  (/  ( U )  UyUyyy  +  ( U )  U  J/J/J/J/ )  |  ( IEj  ,y  J  ) 


such  that 


2(i)  _  2(i) 

J  i+\,i 


Ax 

(7(1)  -  0(1) 


Ay 


/(w)x|(^,%)  +0(Ax5), 
5(«)vl(*i,yJ0  +  0(Ay5), 


and  the  values  of  ufj+±  in  the  monotone  fluxes  /+(•,•)  (approximating  f(u)) 

and  /i2 (• ,  •)  (approximating  g(u))  are  obtained  by  one-dimensional  WENO  interpolation 
in  the  same  procedure  as  in  the  one-dimensional  case.  The  values  of  u  and  all  its 
spatial  derivatives  in  all  the  remaining  terms  are  approximated  by  simple  one-dimensional 
polynomial  interpolations  on  the  same  stencil  as  the  WENO  approximation,  namely  on 
S  =  {xj_ 2,  •  •  • ,  xi+3}  for  fixed  y3  or  on  {yj_ 2, . . . ,  2/7+3}  for  fixed  xt.  This  part  is  also  the 
same  as  in  the  one- dimensional  case. 
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For  the  second  time  derivative,  u"  =  ( f'{u)2ux  +  f'(u)g'(u)uy)x  +  (g'(u)f'(u)ux  + 
g'(u)2uy)y,  we  can  also  approximate  them  in  a  dimension- by-dimension  fashion.  For  the 


fifth  order  scheme,  the  numerical  fluxes  have  the  following  form: 

Ar2 

fl+. ij-  =  -  ( f(u)2ux  +  f'(u)g'(u)uy) \(xi+x,yj)  +  —(2f\u)f"\u)(uxf  +  2  (f"  {u))2{uxf 
+  6f(u)f"(u)uxuxx  +  (. f\u))2uxxx  +  f"\u)g'{u){ux)2uy  +  2  f"{u)g"{u){ux)2uy 
+  f'[u)g"'{u)  {ux)2uy  +  f"(u)g\u)uxxuy  +  f'(u)g"(u)uxxuy 
+  V"{u)g\u)  ^x^xy  +  2f'(u)g"(u)  ^x^xy  +  /W(«)  'U’xxy)  |  (a?  i  ,t/j) 

=  -  ( g\u)f'(u)ux  +  (^(«))2«y)l(*4,yJ+i)  +  ^(V(«)/'(«)K)3  +  2(g" (u))2 (uy)3 
+  Qg'  [u)g"  [u)uyuyy  +  {g' {u))2uyyy  +  g”\u)f\u)ux{uy )2  +  2  g"  (u)  f”  (u)ux{uy)2 
+  g\u)f"\u)ux(uy)2  +  g"{u)f'(u)uxUyy  +  g'(u)  f"  (u)uxuyy 
+  2g'\u)f\u) 

^y^xy  +  2  g'(u)f"(u)  V'y^'xy  "b  fl1  (M)/  (M)'u3:yy)|(a;i,y  i) 


which  satisfies 

?(2)  _  ?(2) 

Ax 


~{f\u)2ux  +  f'{u)g\u)uy)x\{xuyj)  +  0(Ax4)  +  0(Ay4) 


A  (2)  _  A  (2) 

=  -(»'(“)/'(“)“«  +  9'(“)2“»)»lfe.w)  +  0(A**)  +  0(Ay4). 

All  the  spatial  derivatives  can  be  obtained  by  using  simple  one-  or  two-dimensional 
polynomial  interpolations  on  the  same  stencil  as  the  WENO  approximation.  Two- 
dimensional  interpolation  on  the  stencil  S  =  {x*_ 2, . . . ,  Xi+ 3}  x  {yj- 2, . . . ,  ^+3}  is  needed 
to  obtain  approximations  of  the  mixed  derivatives.  This  part  is  different  from  the  one¬ 
dimensional  case.  The  higher  order  time  derivatives  can  be  approximated  in  a  similar 
way. 

Once  the  fluxes  for  all  the  time  derivatives  have  been  obtained,  we  can  form  the  final 
numerical  flux  as  the  sum  of  all  these  fluxes 

1  =  a  1)  ,  ^  m  ,  m  ,  m 

Jl+2  >•?  2  6  24  ■'*+§>.? 

_  a(i)  At  „  (2)  At2  a  (3)  At3  „  (4) 

^4+1  -  ^-+1  +  -y^j+1  + 
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and  obtain  the  final  scheme 


u 


n+1 

hi 


At 


=  u 


hj 


Ax  fi-hi)  A„ 


At 

Ay 


which  is  fifth  order  accurate  in  space  and  fourth  order  accurate  in  time. 

For  the  system  case,  WENO  interpolation  with  local  characteristic  decomposition  is 
still  needed  for  the  computation  of  uf+i  .  and  ufj+ 1  in  the  two-point  fluxes  hi(-,  •)  and 
This  part  is  the  same  as  in  the  one- dimensional  case  and  uses  the  information 
of  the  Jacobians  f'(u )  and  g'{u)  in  the  x  and  y  directions  respectively.  For  all  the  other 
terms,  neither  local  characteristic  decomposition  nor  WENO  interpolation  is  needed. 
We  can  use  the  interpolation  polynomials  as  described  above  for  the  scalar  case  in  each 
component  of  the  solution  u. 


4  Numerical  results 


In  this  section,  we  present  the  results  of  our  numerical  experiments  for  the  schemes 
described  in  the  previous  section,  with  fourth  order  in  temporal  and  fifth  order  in  spatial 
accuracy.  In  some  examples  we  also  compare  with  the  Lax-Wendroff  type  finite  difference 
WENO  scheme  in  [8],  also  with  fourth  order  in  temporal  and  fifth  order  in  spatial 
accuracy. 

Example  1.  We  first  test  the  accuracy  of  the  scheme  on  the  linear  scalar  problem 

ut  +  ux  =  0  (34) 

with  the  initial  condition  u(x,  0)  =  sin(7rx)  and  2-periodic  boundary  condition.  The  exact 
solution  is  u(x,t)  =  sin(7r(x  —  t)).  In  Table  1,  we  show  the  errors  at  time  t  —  2.  We 
can  observe  fifth  order  accuracy  until  quite  refined  meshes,  despite  of  the  formal  fourth 
order  temporal  accuracy,  supporting  our  claim  before  that  spatial  errors  are  usually 
dominating  for  modest  meshes. 
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Table  1:  Accuracy  on  ut  +  ux  =  0  with  u(x,  0)  =  sin(7nr)  at  t  =  2. 


Nx 

L\  errors 

order 

Loo  error 

order 

10 

2.38E-02 

- 

3.67E-02 

- 

20 

1.12E-03 

4.41 

1.99E-03 

4.20 

40 

3.45E-05 

5.02 

6.64E-05 

4.91 

80 

1.07E-06 

5.00 

2.16E-06 

4.94 

160 

3.35E-08 

5.00 

6.51E-08 

5.05 

320 

1.05E-09 

5.00 

1.95E-09 

5.06 

640 

3.25E-11 

5.01 

5.71E-11 

5.10 

Example  2.  We  test  our  scheme  on  the  nonlinear  Burgers  equation 

Ut  +  (y)  =  0  (35) 

with  the  initial  condition  u(x,  0)  =  0.5  +  sin(7rx)  and  a  2-periodic  boundary  condition. 
All  the  three  monotone  fluxes  mentioned  in  Section  3.1,  namely  the  Godunov  flux,  the 
Engquist-Osher  flux  and  the  Lax-Friedrichs  flux,  are  tested. 

When  t  =  0.5/7 r  the  solution  is  still  smooth.  The  errors  for  all  three  monotone  fluxes 
are  shown  in  Tables  2  and  3  respectively.  For  comparison,  we  also  show  the  errors  by 
the  Lax-Wendroff  type  finite  difference  WENO  scheme  in  [8]  with  the  same  fourth  order 
in  temporal  and  fifth  order  in  spatial  accuracy  in  Table  2.  We  observe  that  all  schemes 
achieve  fifth  order  accuracy,  despite  of  the  formal  fourth  order  temporal  accuracy,  once 
again  supporting  our  claim  before  that  spatial  errors  are  usually  dominating  for  modest 
meshes.  For  this  test  case,  there  is  very  little  difference  in  the  errors  corresponding  to 
the  three  different  monotone  fluxes.  We  do  observe  that  the  errors  of  the  WENO  scheme 
in  [8]  are  significantly  bigger  than  those  for  the  schemes  designed  in  this  paper  in  Tables 
2  and  3  when  compared  on  the  same  mesh,  indicating  that  the  narrower  effective  stencil 
and  flexibility  in  choosing  numerical  fluxes  are  helping  to  improve  the  magnitude  of  the 
errors. 

We  plot  the  solution  with  80  grid  points  at  time  t  =  1.5/7T  in  Figure  1,  when  a 
shock  has  already  appeared  in  the  solution.  We  can  see  that  for  all  the  three  monotone 
fluxes,  the  schemes  simulate  the  solution  with  good  resolution  and  without  noticeable 
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oscillations  near  the  shock. 


Table  2:  Accuracy  on  the  Burgers  equation  with  the  Lax- Friedrichs  flux  and  the  WENO 
scheme  in  [8].  u(x,  0)  =  0.5  +  sin(7rx)  at  t  =  0.5/V. 


Lax-Friedrichs  flux 

WENO  scheme  in  [8] 

Nx 

L\  errors 

order 

Loo  error 

order 

L\  errors 

order 

Lqo  error 

order 

10 

4.57E-03 

- 

1.42E-02 

- 

3.65E-02 

- 

6.85E-02 

- 

20 

4.85E-04 

3.24 

2.33E-03 

2.61 

4.44E-03 

3.04 

1.26E-02 

2.45 

40 

2.59E-05 

4.23 

2.38E-04 

3.29 

2.55E-04 

4.12 

1.06E-03 

3.57 

80 

1.37E-06 

4.23 

1.14E-05 

4.38 

9.99E-06 

4.67 

5.11E-05 

4.38 

160 

5.97E-08 

4.52 

9.90E-07 

3.53 

3.76E-07 

4.73 

1.69E-06 

4.92 

320 

1.80E-09 

5.05 

3.69E-08 

4.74 

1.15E-08 

5.03 

7.62E-08 

4.47 

640 

4.07E-11 

5.47 

3.97E-10 

6.54 

2.96E-10 

5.28 

1.46E-09 

5.70 

Table  3:  Accuracy  on  the  Burgers  equation  with  the  Godunov  flux  and  Engquist-Osher 
flux.  u(x,  0)  =  0.5  +  sin(7nr)  at  t  =  0.5/7 r. 


Godunov  flux 

Engquist-Osher  flux 

Nx 

Li  errors 

order 

Lqo  error 

order 

Li  errors 

order 

Lqo  error 

order 

10 

4.59E-03 

- 

1.41E-02 

- 

4.59E-03 

- 

1.41E-02 

- 

20 

4.84E-04 

3.24 

2.32E-03 

2.60 

4.84E-04 

3.24 

2.32E-03 

2.60 

40 

2.59E-05 

4.23 

2.38E-04 

3.29 

2.59E-05 

4.23 

2.38E-04 

3.29 

80 

1.37E-06 

4.23 

1.14E-05 

4.38 

1.37E-06 

4.23 

1.14E-05 

4.38 

160 

5.97E-08 

4.52 

9.90E-07 

3.53 

5.97E-08 

4.52 

9.90E-07 

3.53 

320 

1.80E-09 

5.05 

3.69E-08 

4.74 

1.80E-09 

5.05 

3.69E-08 

4.74 

640 

4.07E-11 

5.47 

3.97E-10 

6.54 

4.07E-11 

5.47 

3.97E-10 

6.54 

Example  3.  We  now  consider  the  one-dimensional  Euler  system  of  compressible  gas 
dynamics: 


&  +  /(£)*  =  0  (36) 

where 

f  =  (p,fm,E)T,  f(£)  =  (pu,pu2 +  p,u(E +  p))T.  (37) 

Here  p  is  the  density,  u  is  the  velocity,  E  is  the  total  energy,  and  p  is  the  pressure, 
which  is  related  to  the  total  energy  by  E  —  +  \pu2  with  7  =  1.4  for  air.  The  initial 

condition  is  set  to  be  p(x,  0)  =  1+0.2  sin(7rx),  u(x,  0)  =  0.7,  p(x,  0)  =  1,  with  a  2-periodic 
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Figure  1:  The  Burgers  equation  with  w(x,0)  =  0.5  +  sin(7nr).  t  =  1.5/7 r.  IV  =  80  grid 
points. 


boundary  condition.  The  exact  solution  is  p(x,t )  =  1  +  0.2sin(7r(a:  —  ut)),  u(x,t)  =  0.7 
and  p(x ,  t)  —  1.  The  errors  and  numerical  orders  of  accuracy  for  the  density  p  are  shown 
in  Table  4.  We  can  see  that  the  schemes  with  with  both  the  Lax-Friedrichs  flux  and 
the  HLLC  flux  can  achieve  the  designed  order  of  accuracy,  while  the  error  magnitude  is 
smaller  for  the  HLLC  flux  than  for  the  Lax-Friedrichs  flux  on  the  same  mesh. 


Table  4:  Errors  and  orders  of  accuracy  for  the  density  p  on  the  one-dimension  Euler 
equation  at  t  =  2.  Lax-Friedrichs  flux  and  HLLC  flux. 


Lax-Friedrichs  flux 

HLLC  flux 

Nx 

L\  errors 

order 

Loo  error 

order 

L\  errors 

order 

Loo  error 

order 

10 

9.20E-03 

- 

1.30E-02 

- 

3.46E-03 

- 

5.58E-03 

- 

20 

4.84E-04 

4.25 

7.74E-04 

4.07 

1.55E-04 

4.48 

2.93E-04 

4.25 

40 

1.53E-05 

4.98 

2.80E-05 

4.79 

4.82E-06 

5.01 

1.00E-05 

4.87 

80 

4.77E-07 

5.00 

8.89E-07 

4.98 

1.50E-07 

5.01 

3.06E-07 

5.04 

160 

1.47E-08 

5.02 

2.67E-08 

5.06 

4.63E-09 

5.02 

8.38E-09 

5.19 

320 

4.42E-10 

5.06 

7.52E-10 

5.15 

1.39E-10 

5.06 

2.37E-10 

5.14 

640 

1.19E-11 

5.21 

1.99E-11 

5.24 

3.72E-12 

5.22 

6.44E-12 

5.20 

Example  4.  Next,  we  solve  the  same  one- dimensional  Euler  equations  with  a  Riemann 
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initial  condition  for  the  Sod’s  problem 


f  (1,0,1),  if  x  <  0 
\  (0.125,0,0.1),  if  x  >  0 


and  the  Lax’s  problem 


f  (0.445,0.698,3.528),  if  x  <  0 
\(0.5,0,0.571),  if  x  >  0. 


(38) 


(39) 


We  plot  the  computed  density  at  t  —  1.3  for  both  problems  with  200  grid  points.  We  can 
see  the  results  of  both  problems  are  well  resolved  and  non-oscillatory  with  Lax-Friedrichs 
flux  and  HLLC  flux. 


(a)  The  Sod’s  problem  (b)  The  Lax’s  problem 

Figure  2:  The  Sod’s  and  Lax’s  problems,  t  =  1.3.  200  grid  points  with  Lax-Friedrichs 
flux  and  HLL  C  flux. 


Example  5.  We  plot  the  numerical  solution  of  the  blast  wave  problem.  We  are  still 
solving  the  Euler  Equations  (36)- (37),  with  the  initial  condition 


P{x,0) 


p(x 

,  OH  1.0, 

0 

<  X 

< 

1 

u(x 

5  0)  =  0.0, 

0 

<  X 

< 

1 

1 

f  1000, 

if 

o  < 

X 

< 

:  0.1 

= 

0.01, 

if 

0.1 

< 

X 

<  0.9 

1 

r - 

o 

o 

if 

0.9 

< 

X 

<  1. 

(40) 
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We  compare  numerical  solution  with  a  reference  solution,  which  is  obtained  with  the 
regular  WENO  scheme  in  [4]  with  16000  grid  points  and  hence  can  be  considered  as  the 
exact  solution.  We  use  both  Lax-Friedrichs  flux  and  HLLC  flux  for  the  problem.  800 
grid  points  are  used.  Our  numerical  results  seem  to  agree  with  the  reference  solution 
very  well. 


(a)  Density 


(b)  Pressure 


(c)  Velocity 

Figure  3:  Blast  wave  problem.  Lines  are  from  the  reference  solution.  Circles  are  the 
numerical  solution  with  Lax-Friedrichs  flux,  and  pluses  are  the  numerical  with  HLLC 
flux:  (a)  density;  (b)  pressure;  (c)  velocity. 


Example  6.  We  now  consider  the  following  two-dimensional  linear  equation 

Ut  +  aux  +  buy  =  0  (41) 
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with  the  initial  condition  u(x,  y,  0)  =  sin(7r(x  +  r/))  and  a  2-periodic  boundary  condition, 
where  a  and  b  are  constants.  Here,  we  take  a  =  1,  b  =  —2,  and  the  final  time  t  =  2.  To 
avoid  possible  error  cancellations  due  to  symmetry,  we  use  different  mesh  sizes  in  the  x 
and  y  directions.  Errors  are  shown  in  Table  5.  We  can  clearly  observe  that  the  scheme 
achieves  the  designed  fifth  order  accuracy. 


Table  5:  Accuracy  on  the  linear  equation  ut  +  ux  —  2uy  =  0  with  u(x,y,  0)  =  sin(7r(x  +  7/)) 
at  t  =  2. 


NX  X  Ny 

L\  errors 

order 

Loo  error 

order 

8  x  12 

4.97E-02 

- 

7.23E-02 

- 

16  x  24 

5.70E-03 

3.12 

1.02E-02 

2.83 

32  x  48 

1.70E-04 

5.07 

3.42E-04 

4.90 

64  x  96 

4.72E-06 

5.17 

1.01E-05 

5.08 

128  x  192 

1.41E-07 

5.07 

3.07E-07 

5.04 

256  x  384 

4.32E-09 

5.02 

9.16E-09 

5.06 

Example  7.  Next,  we  solve  the  two-dimensional  nonlinear  Burgers  equation 


with  the  initial  condition  u(x,  0)  =  0.5  +  sin(7r(x  +  y)/ 2)  and  a  4-periodic  boundary 
condition.  When  t  =  0.5/7T  the  solution  is  still  smooth,  and  we  give  the  errors  and  order 
of  accuracy  in  Table  6  with  the  Lax-Friedrichs  flux.  Clearly,  the  scheme  delivers  the 
designed  order  of  accuracy.  For  comparison,  we  also  show  the  errors  by  the  Lax-Wendroff 
type  finite  difference  WENO  scheme  in  [8]  with  the  same  fourth  order  in  temporal  and 
fifth  order  in  spatial  accuracy  in  Table  6.  We  again  observe  that  the  errors  of  the  WENO 
scheme  in  [8]  are  significantly  bigger  than  those  for  the  schemes  designed  in  this  paper  in 
Table  6  when  compared  on  the  same  mesh,  indicating  that  the  narrower  effective  stencil 
is  helping  in  improving  the  magnitude  of  the  errors. 

We  plot  the  results  at  t  =  1.5/7 r  when  shocks  have  already  appeared  in  Figure  4  with 
160  x  160  grid  points.  On  the  left  we  are  showing  a  one-dimensional  cut  at  x  —  y  of  the 
solution,  and  on  the  right  we  are  showing  the  surface  of  the  computed  solution.  Clearly, 
the  shocks  are  well  resolved  without  spurious  oscillations. 
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2  2 

Table  6:  Accuracy  on  the  Burgers  equation  ut  +  (\)x  +  {\)y  —  0  with  Lax- Friedrichs 
flux  and  the  WENO  scheme  in  [8].  u(x,  y,  0)  =  0.5  +  sin(7r(x  +  y)/2)  at  t  —  0.5/n. 


Lax-Friedrichs  flux 

WENO  scheme  in  [8] 

NX  X  Ny 

L\  errors 

order 

Loo  error 

order 

L\  errors 

order 

Loo  error 

order 

8  x  12 

7.51E-03 

- 

2.35E-02 

- 

2.23E-02 

- 

6.44E-02 

- 

16  x  24 

9.86E-04 

2.93 

6.96E-03 

1.75 

2.98E-03 

2.90 

1.45E-02 

2.15 

32  x  48 

8.23E-05 

3.58 

7.30E-04 

3.25 

2.90E-04 

3.36 

1.65E-03 

3.13 

64  x  96 

4.26E-06 

4.27 

4.09E-05 

4.16 

8.58E-06 

5.08 

8.08E-05 

4.36 

128  x  192 

1.84E-07 

4.53 

1.62E-06 

4.65 

3.27E-07 

4.71 

2.57E-06 

4.97 

256  x  384 

6.24E-09 

4.88 

7.20E-08 

4.49 

1.01E-08 

5.01 

1.33E-07 

4.28 

Figure  4:  The  two-dimensional  Burgers  equation.  u(x,Q)  =  0.5  +  sin(7r(x  +  y)/ 2).  t  = 
1.5/7 r.  160  x  160  grid  points.  Left:  a  cut  of  the  solution  at  x  =  y,  where  the  solid  line 
is  the  exact  solution  and  the  circles  are  the  computed  solution;  right:  the  surface  of  the 
solution. 


Example  8.  We  now  consider  two-dimensional  Euler  systems  of  compressible  gas  dy¬ 
namics: 


+  /(£)*  +  9(Oy  —  0 


(43) 
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with 


f  =  ( P,pu,pv,E)T , 

/(O  =  (pit,  pw2  +  P,  puv,  u(E  +  p))r, 
s(f )  =  (/w,  pun,  pv2  +  p,  v(E  +  p))T. 

Here,  p  is  the  density,  (u,  v )  is  velocity,  E  is  the  total  energy,  and  p  is  the  pressure, 
related  to  the  total  energy  E  by  E  =  +  \p(u2  +  n2),  with  7  =  1.4.  The  initial 

condition  is  set  to  be  p(x,  y,  0)  =  1  +  0.2  sin  (77  (x  +  y)),  u(x,  y ,  0)  =  0.7,  v(x,  y,  0)  =  0.3, 
p(x,y,  0)  =  1,  with  a  2-periodic  boundary  condition.  The  exact  solution  is  p(x,y,t )  = 
1  +  0.2  sin(7r(x  +  y  —  (u  +  v)t)),  u  =  0.7,  v  =  0.3,  p  =  1.  The  errors  and  order  of  accuracy 
for  the  density  are  shown  in  Table  7.  We  can  see  that  the  scheme  achieves  the  designed 
order  of  accuracy. 


Table  7:  Accuracy  on  the  two-dimensional  Euler  equation  at  t  =  2.  Errors  of  the  density 

P- 


NX  X  Ny 

L\  errors 

order 

Loo  error 

order 

8  x  12 

2.29E-02 

- 

3.28E-02 

- 

16  x  24 

1.47E-03 

3.97 

2.22E-03 

3.89 

32  x  48 

5.15E-05 

4.83 

9.03E-05 

4.62 

64  x  96 

1.61E-06 

5.00 

2.99E-06 

4.92 

128  x  192 

4.98E-08 

5.01 

9.20E-08 

5.02 

256  x  384 

1.51E-09 

5.04 

2.58E-09 

5.15 

Example  9.  Next,  we  use  the  following  problem  to  illustrate  further  the  high  order 
accuracy  of  the  method.  Consider  the  following  idealized  problem  for  the  Euler  equation 
in  2D:  The  mean  flow  is  p  —  1,  p  —  1,  and  (u,v)  =  (1, 1)  (diagonal  flow).  We  add,  to 
this  mean  flow,  an  isentropic  vortex  (perturbation  in  (u,v)  and  the  temperature  T  =  2, 
no  perturbation  in  the  entropy  S  —  ^7): 

(Su,Sv)  =^-e0-5{1~r2\-y,x) 

Z7T 

ST  =  -  (7  ~  y  e1-’ 


877 T2 


5S  =0 


where  (cc,y)  =  (x  —  5 ,y  —  5),  r2  =  x2  +  y2,  and  the  vortex  strength  e  =  5  [11].  Since 
the  mean  flow  is  in  the  diagonal  direction,  the  vortex  movement  is  not  aligned  with  the 
mash  direction.  It  is  clear  the  the  exact  solutions  is  just  the  passive  convection  of  the 
vortex  with  the  mean  velocity.  The  computational  domain  is  taken  as  [0, 10]  x  [0,10], 
extended  periodically  in  both  directions.  The  simulation  is  performed  until  t  =  0.2.  The 
errors  and  orders  of  accuracy  are  shown  in  Table  8. 


Table  8:  Accuracy  on  the  Vortex  evolution  at  t  =  0.2.  Errors  of  the  density  p. 


NX  X  Ny 

L  i  errors 

order 

Loo  error 

order 

40  x  40 

5.52E-03 

7.50 

1.98E-03 

8.98 

80  x  80 

2.90E-04 

4.25 

1.41E-04 

3.81 

160  x  160 

6.84E-06 

5.41 

2.57E-06 

5.78 

320  x  320 

1.96E-07 

5.13 

5.53E-08 

5.54 

Example  10.  Finally,  we  consider  the  double  Mach  reflection  problem.  The  computa¬ 
tional  domain  for  this  problem  is  chosen  to  be  [0,4]  x  [0, 1].  The  reflection  wall  lies  at 
the  bottom  of  the  computational  domain  starting  from  x  —  |.  Initially  a  right- moving 
Mach  10  shock  is  positioned  at  x  —  |,  y  —  0  and  makes  a  60°  angle  with  the  x-axis.  For 
the  bottom  boundary,  the  exact  post-shock  condition  is  imposed  for  the  part  from  x  =  0 
to  x  =  |  and  a  reflective  boundary  condition  is  used  for  the  rest.  At  the  top  boundary 
of  the  computational  domain,  the  flow  values  are  set  to  describe  the  exact  motion  of  the 
Mach  10  shock.  The  initial  pre-shock  condition  is 

(p,p,u,v)  =  ( 8,  116.5,  8.25 cos(30°),  -8.25 sin(30°))  (44) 

and  the  post-shock  condition  is 

(p,P,u,v)  =  (1.4,  1,  0,  0)  (45) 

We  compute  the  solution  up  to  t  =  0.2.  Uniform  meshes  with  grid  points  960  x  239  are 
used.  In  Figure  5  we  show  30  equally  spaced  density  contours  from  1.5  to  22.7,  and  a 
“zoomed-in”  graph  is  provided  in  Figure  6.  We  can  see  that  our  results  agree  well  with 
the  results  in  [8]. 
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Figure  5:  Double  Mach  reflection  problem.  30  equally  spaced  contours  form  1.5  for  22.1 
for  the  density  p. 


Figure  6:  Zoomed-in  figure.  Double  Mach  reflection  problem.  30  equally  spaced  contours 
form  1.5  for  22.7  for  the  density  p. 

5  Concluding  remarks 

We  have  discussed  in  depth  an  alternative  formulation  for  constructing  high  order  con¬ 
servative  finite  difference  schemes.  In  this  formulation,  the  high  order  interpolation 
procedure  is  applied  to  the  solution  itself  rather  than  to  the  flux  functions.  Even  though 
the  algebra  and  algorithm  complexity  are  increased  compared  with  the  standard  finite 
difference  formulation,  this  alternative  formulation  does  have  several  advantages,  includ¬ 
ing  the  flexibility  in  using  any  monotone  fluxes  in  the  scalar  case  and  any  approximate 
Riemann  solvers  in  the  system  case,  and  the  easiness  in  applying  the  Lax-Wendroff  time 
discretization,  resulting  in  a  scheme  with  a  narrow  stencil  which  is  the  same  as  that  for 
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the  semi-discrete  scheme.  We  demonstrate  these  two  advantages,  and  discuss  the  second 
advantage  in  depth.  We  describe  the  Lax-Wendroff  time  discretization  of  this  alternative 
formulation  of  high  order  conservative  finite  difference  schemes  in  detail,  using  fifth  order 
WENO  interpolation  for  the  leading  term  and  fixed  stencil  polynomial  interpolation  for 
the  remaining  terms.  Temporal  accuracy  is  kept  at  fourth  order  in  the  extensive  nu¬ 
merical  results,  including  one-  and  two-dimensional  scalar  problems  and  systems.  These 
numerical  examples  verify  the  performance  of  the  designed  schemes  in  their  designed 
high  order  of  accuracy  and  non-oscillatory  shock  transitions  for  discontinuous  solutions. 
Another  advantage  of  this  alternative  formulation,  namely  the  easier  maintenance  of 
free-stream  solutions  in  multi-dimensional  curvilinear  coordinates,  will  be  explored  in 
future  work. 
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